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1.  Introduction 


An  important  extension  of  the  original  constant-temperature  dissipative  particle  dynamics  (DPD) 
method  (1 ,  2)  that  imposes  constant  energy  (DPD-E)  conditions  was  developed  by  Bonet  Avalos 
and  Mackie  (3)  and  later  independently  by  Espanol  ( 4 ).  The  DPD-E  method  includes  an 
additional  equation-of-motion  that  provides  a  dynamic  depiction  of  the  internal  state  of  a  coarse- 
grain  particle.  Consequently,  particles  are  allowed  to  exchange  both  momentum  and  heat,  thus 
satisfying  total  energy  and  momentum  conservation. 

Numerical  integration  of  the  equations  of  motion  (EOMs)  is  a  key  consideration  when  applying 
the  DPD  method  because  the  stochastic  component  requires  special  attention.  Efficient  and 
accurate  integration  schemes  are  required  to  maintain  reasonable  time  step  sizes,  thus  allowing 
for  the  simulation  of  mesoscale  events.  Moreover,  the  advent  of  conservative  forces  extending 
beyond  purely  repulsive  models  that  contain  attractive  character  further  supports  the  need  for 
effective  integration  schemes.  However,  the  integration  is  a  nontrivial  task  due  to  the  pairwise 
coupling  of  particles  through  the  random  and  dissipative  forces  (5).  For  example,  the  most 
widely  used  modified  velocity- Verlet  algorithm  ( 6)  works  reasonably  well  for  the  constant- 
temperature  DPD  method,  but  for  the  DPD-E  method  it  requires  time-step  values  that  are  several 
orders  of  magnitude  smaller  than  for  constant-temperature  DPD  (7,  8). 

Furthermore,  self-consistent  solutions  are  often  necessary  because  the  dissipative  forces  and  the 
relative  velocities  of  the  particles  are  interdependent,  where  the  considerable  computational  cost 
associated  with  this  task  has  motivated  the  development  of  various  integration  schemes.  Recent 
independent  studies  by  Nikunen  et  al.  (5)  and  Chaudhri  and  Lukes  (9)  assessed  the  quality  and 
performance  of  several  applicable  integration  schemes  for  constant-temperature  DPD,  where  the 
Shardlow-splitting  algorithm  (SSA)  (10)  was  identified  as  the  best-performing  approach.  In 
recent  work  by  our  group  (11, 12),  a  comprehensive  description  of  numerical  integration 
schemes  based  upon  the  SSA  was  presented  for  both  the  isothermal  and  isothermal,  isobaric 
DPD  methods.  The  original  SSA  formulated  for  systems  of  equal-mass  particles  was  extended  to 
systems  of  unequal-mass  particles.  Both  a  velocity- Verlet  scheme  and  an  implicit  scheme  were 
formulated  for  the  integration  of  the  fluctuation-dissipation  contribution,  where  the  velocity- 
Verlet  scheme  consistently  performed  better. 

The  SSA  decomposes  the  EOM  into  differential  equations  that  correspond  to  the  deterministic 
dynamics  and  elementary  stochastic  differential  equations  that  correspond  to  the  stochastic 
dynamics.  In  the  original  SSA  formulation,  both  types  of  differential  equations  are  integrated  via 
the  velocity- Verlet  algorithm  (10),  where  the  stochastic  dynamics  are  additionally  solved  in  an 
implicit  manner  that  conserves  linear  momentum.  Previously,  the  SSA  had  only  been  applied  to 
the  isothermal  (10, 11)  and  isothennal,  isobaric  (12)  DPD  methods. 
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In  this  work,  we  formulate  the  SSA  for  the  DPD-E  method,  where  we  verify  the  DPD-E  variant 
using  both  the  standard  DPD  fluid  (pure  and  binary  mixtures)  ( 6)  and  a  coarse-grain  solid  model 
(13).  For  both  models,  we  further  verify  the  DPD-E  variant  by  instantaneously  heating  a  slab  of 
particles  in  the  simulation  cell.  We  monitor  the  evolution  of  the  temperature  and  pressure  as  the 
system  approaches  an  equilibrated  state  while  maintaining  constant  energy.  For  completeness, 
the  derivations  of  the  Fokker-Planck  equation  (FPE)  and  the  fluctuation-dissipation  theorem 
(FDT)  are  included  in  appendix  A. 


2.  Formulations  of  DPD  at  Fixed  Total  Energy  Using  Shardlow-Like 
Splitting  Numerical  Discretization 


2.1  General  Formulation  of  DPD 

DPD  particles  are  defined  by  a  mass  mi ,  position  r, ,  and  momentum  p,  .  The  particles  interact 
with  each  other  via  a  pairwise  force  Fy-  that  is  written  as  the  sum  of  a  conservative  force  Ff , 
dissipative  force  V'5 ,  and  random  force  Ff  : 


=  f,c+f"+f;. 


(i) 


Ff  is  given  as  the  negative  derivative  of  a  coarse-grain  potential,  u(i;  G ,  expressed  as 


du 


CG 


r„ 


dru  rv 


(2) 


where  r  ;  =  r  -  r  is  the  separation  vector  between  particle  i  and  particle  j ,  and  rtj  =  r,y  .  The 

remaining  two  forces,  F^  and  Fy  ,  can  be  interpreted  as  a  means  to  compensate  for  the  degrees 

of  freedom  neglected  by  coarse  graining.  The  conservative  force  is  not  specified  by  the  DPD 
formulation  and  can  be  chosen  to  include  any  forces  that  are  appropriate  for  a  given  application, 
including  multibody  interactions  (e.g.,  [13-15]).  F,.f  and  F*  are  defined  as 


(3) 


and 


TJ_ 
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where  between  particle  i  and  j ,  and  cry  are  the  friction  coefficient  and  noise  amplitude  , 
respectively,  v ,  =  —  -  —  and  W..  are  independent  Wiener  processes  such  that  Wu  =  W n  .  The 


m,  m , 


J 


weight  functions  coD{r)  and  t oR{r )  vanish  for  r>rc ,  where  r  is  the  cutoff  radius. 


Note  that  Fc  is  completely  independent  of  F/J  and  Ff  ,  while  F"  and  FA'  are  not  independent 

but  rather  are  coupled  through  a  fluctuation-dissipation  relation.  This  coupling  arises  from  the 
requirement  that  in  the  thermodynamic  limit,  the  system  samples  the  corresponding  probability 
distribution.  The  necessary  conditions  can  be  derived  using  an  FPE;  these  conditions  are 
presented  in  appendix  A  for  the  DPD-E  variant. 


2.2  Constant-Energy  DPD 

The  constant-energy  DPD  approach  was  developed  by  Bonet  Avalos  and  Mackie  (3),  and  shortly 
thereafter  also  by  Espanol  ( 4 ).  To  conserve  energy  in  a  DPD  simulation,  an  additional  variable  is 
introduced  that  characterizes  the  internal  state  of  the  particles.  An  internal  energy  ut  (restricted 

to  ui  >  0  always)  is  associated  with  each  DPD  particle,  which  accounts  for  the  energy  absorbed 

or  released  by  the  atomic  internal  degrees  of  freedom  that  have  been  coarse-grained  into  the 
DPD  particle.  The  total  energy  of  the  system  is  conserved  since  the  kinetic  energy  lost/gained  by 
the  dissipative  and  random  interactions  is  absorbed/released  by  this  particle  internal  energy. 
Taken  along  with  uj  is  an  associated  mesoscopic  entropy  y  =  s(uj),  which  is  a  monotonously 

1  c)s  ■ 

increasing  function  of  ui ,  so  that  the  internal  temperature  0:  is  defined  as  —  =  — >0  (16).  A 

0i  dui 

mesoparticle  equation  of  state  relating  0I  and  ui  is  therefore  required.  It  is  convenient  to  write 
the  variation  of  d ui  as  the  sum  of  two  contributions  that  correspond  to  the  mechanical  work  done 
on  the  system  d u™ech  and  the  heat  conduction  between  particles  d u-ond  ,  i.e.,  as 
dw,  =  d u™ech  +  d uc°"d  .  The  dynamics  of  the  system  is  then  governed  by  the  following  equations- 
of-motion  (EOMs): 


dr;  =  —  dt 
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where  coDq(r)  and  coRq{r)  are  weight  functions  vanishing  for  r  >  rc ,  k{j  and  cr  are  the 
mesoscopic  thermal  conductivity  and  the  noise  amplitude  between  particle  i  and  particle  j , 
respectively,  and  dWq  =  -d  Wq  are  the  increments  of  the  Wiener  processes  associated  with 

thennal  conduction.  Equation  5  is  a  generalization  of  the  modified  DPD-E  approach  ( 16)  for 
particles  with  unequal  masses  (17).  In  the  last  expression  of  equation  5,  the  first  term  on  the 
right-hand  side,  which  specifies  the  dissipative  heat  conduction,  can  alternatively  be  expressed  in 
tenns  of  the  difference  of  particle  temperatures  (18). 


Bonet  Avalos  and  Mackie  (5)  demonstrated  that  thermodynamic  consistency  requires  the 
following  fluctuation-dissipation  relations  to  be  satisfied: 

=  2Yijk*®j 
OJD(r)=  [ft/fr)]2 

a]  =  2  kBKy 
0)Dq(r)=[o)Rq(r)f 


(6) 


where  the  relevant  temperature  is 


f 


and  coDq(r)  and  coRq(r)  can  be  chosen 
similar  to  coD{r )  and  t oR{r ),  respectively  (11).  (Bonet  Avalos  and  Mackie  [J]  did  not  prove  that 


the  DPD-E  EOMs  sample  the  microcanonical  ensemble,  rather  they  proved  that  the  relations  in 
equation  6  are  required  for  the  EOMs  to  sample  the  canonical  ensemble,  where  these  relations  are 
independent  of  T  .)  An  outline  of  the  derivation  of  these  FDTs  along  with  the  FPE  is  given  in 
appendix  A.  By  design,  these  EOMs  conserve  total  momentum  P  and  the  total  energy 
E  —  U  +  K  +  y  u, .  Note  that  the  random- force  noise  amplitude  cr  depends  on  the  particle 


internal  temperatures  and  not  on  the  system  temperature  T  as  it  does  in  the  constant-temperature 


DPD  method.  (The  system  temperature  in  a  DPD-E  simulation  is  defined  as 


1  yP,  P/ 

3  NkB  i  mi 


•) 


2.2.1  Numerical  Discretization 

As  part  of  the  founding  work,  Mackie  et  al.  (16)  developed  an  explicit  integration  algorithm  for 
the  DPD-E  approach.  However,  because  of  round-off  error  resulting  from  the  loss  of  mechanical 
energy  during  integration  of  the  work  done  by  the  dissipative  and  random  forces,  the  algorithm 
requires  a  rather  small  At  to  satisfactorily  conserve  the  total  energy.  This  situation  can  be 
improved  by  extending  the  splitting  strategy  developed  by  Stoltz  (17)  to  equation  5.  Stoltz 
introduced  a  set  of  EOMs  that  closely  resembles  DPD-E  except  (1)  they  neglect  thermal 
conduction,  and  (2)  they  do  not  project  the  dissipative  and  random  forces  along  the  separation 
vectors  between  the  particles. 
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In  the  following  paragraphs,  we  present  an  extension  of  Stoltz’s  splitting  algorithm  to  DPD-E. 

As  done  previously  (11, 12),  we  start  by  decomposing  the  EOMs  given  in  equation  5  into 
deterministic  differential  equations  and  elementary  stochastic  differential  equations  (SDEs) 
based  upon  the  conservative  and  fluctuation-dissipation  contributions,  respectively.  The 
conservative  terms  are  identical  to  the  constant-temperature  DPD  formulation, 

dr=  —  At  (7a) 

mt 


and 


dpi=EF(f*. 


J*l 


while  the  fluctuation-dissipation  terms  can  be  expressed  as 
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where  the  superscript  i  -  j  indicates  that  the  variation  of  momenta  is  considered  for  a  pair  of 
interacting  particles  i  and  j  only,  and  AW0  =  d  Wfl  are  the  increments  of  the  Wiener  processes. 

Equation  8b  directly  follows  from  the  introduction  of  the  mechanical  contribution  to  the  internal 
energies  (5,  4),  where  Mackie  et  al.  (16)  showed  that  the  total  energy  of  a  pair  of  interacting 
particles  i  and  j  is  conserved  when 
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Using  equation  8a  along  with  the  definitions  of  F^  and  F*  (equations  3  and  4,  respectively),  we 
can  rewrite  equation  9a  as 
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Equation  9b  together  with  — - -  — - 


mech,i—j 


—  then  leads  to  equation  8b.  As  done  previously 

dt  dt 

(11, 12),  the  stochastic  flow  map  (f)M  can  be  approximated  from  the  first-order  splitting  algorithm 
given  by 


/  idiss  idiss  idiss  idiss  idiss  iC  / 1  a\ 

YAt  =  r  At\\,2  °  YAt:L3  °  ***  °  YAtAJ  °  •••  °  YAt:N-2,N  °  YAt:N-LN  °  Y A f>  (10) 


"At;N-2,N  YAt;N-l,N 
idiss 


where  o  denotes  a  composition  of  operators.  For  each  .  term  at  fixed  internal  temperatures 
(  0I  and  0 . ),  momenta  are  updated  through  the  Shardlow-splitting  algorithm-velocity  Verlet 
(SSA-VV)  approach  based  upon  the  constant-temperature  DPD  formulation  previously  given  (11). 
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where  the  superscript  i  -  j  has  been  omitted  for  notational  simplicity,  and  gtj  =  gjt  is  a  Gaussian 
random  number  with  zero  mean  and  unit  variance  that  is  chosen  independently  for  each  pair  of 

interacting  particles,  and  ju  =  —  +  —  . 

mi  mj 


The  conductive  contribution  to  the  internal  energies  is  updated  using  a  Euler  algorithm.  For  the 
SSA-VV,  this  results  in  the  following  expressions: 
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In  equation  12a,  gq  =  -gt  is  a  Gaussian  random  number  with  zero  mean  and  unit  variance, 

chosen  independently  for  each  pair  of  interacting  particles.  Next,  the  mechanical  contribution  to 
the  internal  energies  is  updated  using 
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In  equations  12a  and  12b,  the  superscript  i  —  j  is  again  omitted  for  notational  simplicity.  The 
total  system  energy  is  exactly  conserved  via  equation  12b.  Finally,  <f)(y  can  be  approximated  by 
the  velocity- Verlet  algorithm  previously  given  (11). 
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The  following  outline  summarizes  a  practical  implementation  of  the  SSA-VV  for  the  DPD-E 
variant. 
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2.  Deterministic  Integration  #1  for  i  =  1,...,  /V 


7?  ru  Va7 

-<V°  ?// 


r?  2 


l  c 

(i)  P/^-Pz  +  yF, 

(ii)  r;  <—  r.  +  Ai — 

//7, 

3.  Conservative  Force  Calculation:  |f,c  j.^ 

4.  Deterministic  Integration  #2  for  i  =  1,...,  /V 
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3.  Computational  Details 


The  SSA-VV  for  the  DPD-E  variant  was  tested  using  both  the  standard  DPD  fluid  ( 6)  and  a 
coarse-grain  solid  model  (13),  where  complete  details  of  the  conservative  forces  for  these  models 
are  given  in  appendix  B.  Both  a  pure  component  case  and  an  equimolar  binary  mixture  were 
tested  for  the  DPD  fluid  model.  System  sizes  for  the  DPD  fluids  and  coarse-grain  solid  were, 
respectively,  N  =  10125  and  13500.  For  these  simulations,  the  following  reduced  units  were 
used:  rc  and  r0  are  the  unit  of  length  for  the  DPD  fluid  and  coarse-grain  solid,  respectively;  the 

mass  of  a  DPD  particle  is  the  unit  of  mass;  and  the  unit  of  energy  is  kBTini ,  where  Tini  is  the 

initial  system  temperature.  Using  these  reduced  units,  we  set  the  maximum  repulsion  between 
particles  i  and  j  as  a,  =  25  for  the  pure  DPD  fluid,  and  as  a,  =  25  and  28  for  the  like  and 

unlike  i  -  j  interactions,  respectively,  for  the  binary  DPD  fluid.  Further,  for  all  cases,  we  set  the 
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noise  amplitude  cr  =  3  .  Next,  we  assume  that  the  internal  degrees  of  freedom  are  purely 
hannonic  and  express  the  coarse-grain  particle  equation  of  state  as  uj  =  Cv  i9i ,  with  heat  capacity 
CV  i  /  kE  =  Cv  /  kB=  60  and  48  for  the  pure  DPD  fluid  and  coarse-grain  solid,  respectively.  Note 

that  these  values  correspond  to  coarse-graining  approximately  20  atoms  and  16  atoms  into  a 
DPD  particle,  respectively  (19).  For  the  binary  DPD  fluid,  we  set  m2  =  10 ml  and  Cv  2  =  10CF , , 

where  CV1  /  kB  =  60 .  Finally,  following  Ripoll  et  al.  (20),  the  mesoscopic  thermal  conductivity 
Ky  is  chosen  as 

<14> 

where  kq  is  the  parameter  controlling  the  thermal  conductivity  of  the  DPD  particles,  which  are 
chosen  as  k0  =  2.80  •  10~4  for  the  pure  DPD  fluid  and  k0  =  1.52  •  10  4  for  the  coarse-grain  solid. 
For  the  binary  fluid,  we  set  sr0  n  =  2.80  •  1 0  4 ,  kq  22  =  2.80  ■  1 0  6 ,  and  k0  12  =  ^/f0J1/f0  22  .  When 
6i  =6 .  =  1 ,  these  particular  values  of  kq  and  Cv  give  /v  =  1 . 


4.  Results 


This  section  is  organized  as  follows.  The  validity  of  the  SSA  integration  algorithm  for  the  DPD- 
E  variant  is  verified  by  considering  equilibrium  and  nonequilibrium  scenarios,  where  results  are 
given  in  subsections  4.1  and  4.2,  respectively.  Finally,  we  briefly  review  the  energy  drift 
associated  with  finite  integration  methods  and  propose  a  simple  strategy  to  minimize  these  drifts 
in  DPD-E  simulations. 

4.1  Test  Case  1:  Equivalence  of  DPD  Variants 

As  a  first  test  of  the  SSA-VV  formulation  for  the  DPD-E  variant,  we  verify  that  it  converges  to 
the  same  equilibrium  properties  when  at  the  same  thermodynamic  conditions  as  a  constant- 
temperature  DPD  simulation. 

4.1.1  DPD  Fluid 

The  benchmark  systems  for  both  the  pure  and  binary  DPD  fluid  cases  are  taken  from  a  constant- 
temperature  DPD  simulation  performed  at  p  =  3  and  T  =  1 ,  and  run  for  tnm  =  3000  and 

At  =  0.03  .  The  following  quantities  were  evaluated  and  are  listed  in  tables  1  (pure  fluid)  and  2 
(equimolar  binary  fluid):  virial  pressure  (Pvir) ,  configurational  energy  per  particle  (it) ,  kinetic 

temperature  (Tkjn) ,  and  self-diffusion  coefficients  D  using  the  Einstein  relation  (21). 
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To  validate  the  constant-energy  SSA-VV  formulation,  DPD-E  simulations  were  performed  at 
conditions  taken  from  the  constant-temperature  DPD  simulation,  i.e.,  V  =  E  =  3375  (L  is  the 
box  length).  The  final  configuration  of  the  constant-temperature  DPD  simulation  is  used  to 
determine  the  imposed  values  of  E  ,  thus  it  is  also  used  as  the  starting  configuration  for  the 
DPD-E  simulation.  For  both  the  pure  and  equimolar  binary  fluid  cases,  the  values  of  ui  were 

initialized  by  setting  ui  =  Cv  Tini ,  where  Tini  =T  =  1,  and  were  carried  out  for  tnm  =  3000  and 

At  -  0.01 .  Analogous  to  microcanonical  MD  simulations,  the  use  of  a  smaller  At ,  with  respect 
to  constant-temperature  DPD  simulations,  is  required  for  proper  conservation  of  E  .  Using 
At  =  0.01 ,  we  observed  a  relative  drift  in  E  no  greater  than  1  •  10  4 .  (For  the  DPD  fluid 
simulations,  reported  relative  drifts  refer  to  an  average  of  relative  drifts  over  time  periods  of 
1000.) 


Comparing  the  DPD-E  results  with  the  constant-temperature  DPD  results  in  tables  1  and  2,  we 
find  excellent  overall  agreement.  For  the  DPD-E  simulations,  the  internal  temperature, 

/  1  N  1  \~* 

Tint)  =  i  — ^ —  )  was  also  evaluated,  where  the  values  of  (Tkm)  and  (Tint)  agree  within 

\N  i=i  °i 


statistical  uncertainties.  (Since  0t  is  defined  as  a  ratio  of  si  and  ut ,  \Tmt)  is  estimated  through  a 

hannonic  average  rather  than  an  arithmetic  average  [16,  77].)  For  the  pure  DPD  fluid,  these 
values  are  approximately  1.5%  lower  than  Tini  =  1 ;  however,  this  discrepancy  is  due  to  the 

fundamental  differences  between  the  constant- temperature  DPD  and  DPD-E  methods. 
Effectively,  the  two  systems  are  different  since  the  imposed  temperature  for  the  constant- 
temperature  DPD  system  should  be  equivalent  to  (Tkin ) ,  while  in  the  DPD-E  system  the  total 

energy  initially  given  to  the  system  is  dynamically  partitioned  among  the  kinetic  and  internal 
energies,  yielding  a  variation  in  the  equilibrium  temperature  with  respect  to  Tmj  (16).  This 

difference  is  of  0(kB  /  Cv )  as  compared  with  unity,  while  an  additional  contribution  of  the  same 
order  arises  from  the  “extra  degree-of-freedom”  due  to  the  fluctuations  in  uj ,  since  the  relevant 


macroscopic  temperature  is  related  to 


.  Hence,  for  a  pure  DPD  fluid  up  to  first  order  in 


kJCv, 


Tkin)  =  (T’nt)  -  Tim 0  / Cv)=  0.983  (16),  in  agreement  with  the  simulated  values  of 


0.985  ±  0.003 .  For  the  equimolar  binary  fluid,  the  simulated  values  are  also  in  agreement  with 
the  estimate  (Tkin )  =  (Tmt  'j  =  0.994  .  The  estimated  value  is  closer  to  Tini  =  1  because  of  the  larger 

value  of  CV  2  that  reduces  the  (kB  /  Cv )  contribution.  Also  note  that  because  of  these  lower  values 
of  (Tkin)  and  (7’nt) ,  the  values  of  (Pvir)  in  tables  1  and  2  slightly  differ  from  (Pvir)  for  constant- 
temperature  DPD. 
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Table  1 .  The  configurational  energy  per  particle  (ll  ^  ,  the  kinetic 
temperature  (Tkjn  'j  ,  the  internal  temperature  (Tjnt  \  ,  the 
virial  pressure  (  P.ir  ^ ,  and  the  self-diffusion  coefficient  D 
determined  from  test  case  1  simulations  of  the  pure  DPD 
fluid.  (?)  denotes  an  ensemble  average,  where  numbers  in 

parentheses  are  uncertainties  calculated  from  block 
averages. 


Variant 

M 

(T«.) 

(t.  ) 

y-  int  / 

(p*) 

D 

DPD 

P  =  3 

4.56(1) 

1.005(8) 

— 

23.65(8) 

0.295(13) 

DPD-E 

P  =  3 

4.54(1) 

0.985(8) 

0.985(3) 

23.61(11) 

0.293(6) 

Table  2.  The  configurational  energy  per  particle  (u') ,  the  kinetic  temperature 

(Tkin  'j  ,  the  internal  temperature  (Tm  ^  ,  the  virial  pressure  (Pvjr  ^  ,  and  the 
self-diffusion  coefficient  D  determined  from  test  case  1  simulations  of 
the  equimolar  binary  DPD  fluid.  (?j  denotes  an  ensemble  average,  where 
numbers  in  parentheses  are  uncertainties  calculated  from  block  averages. 


Variant 

(u) 

(T*) 

(t.  ) 

y  int  / 

(p*) 

A 

A 

DPD 

P  =  3 

4.76(1) 

1.005(8) 

— 

24.79(13) 

0.177(13) 

0.165(13) 

DPD-E 

P  =  3 

4.75(1) 

0.993(8) 

0.994(3) 

24.77(21) 

0.174(9) 

0.161(10) 

Reproducing  equilibrium  averages  is  necessary  but  not  sufficient  proof  that  the  integration 
scheme  is  behaving  properly.  Hence,  as  a  further  demonstration  of  the  quality  of  the  SSA-VV 
and  the  proper  choice  of  At ,  for  the  pure  DPD  fluid,  we  calculated  probability  distributions  for 
Pj ,  ui ,  and  V  for  constant-temperature  DPD  with  At  =  0.03  and  for  DPD-E  with  At  =  0.01 .  We 
compared  the  probability  distribution  for  pi  with  the  corresponding  Maxwell-Boltzmann 
distribution  (21),  while  the  probability  distributions  for  ui  and  V  were  compared  with  those 
obtained  with  a  very  small  At  =  0.001 ,  which  is  more  than  an  order  of  magnitude  smaller  than 
typical  values  of  At  used  and  thus  is  approximated  as  the  “exact”  result.  For  a  special  case  of 
DPD-E  in  the  absence  of  conservative  forces,  an  analytical  fonn  of  the  probability  distribution 
for  il  was  derived  by  Mackie  et  al.  (16)  under  constant-temperature  conditions.  Probability 

distributions  from  constant-temperature  DPD  and  DPD-E  (not  shown  here)  are  in  extremely 
good  agreement. 
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4.1.2  Coarse-Grain  Solid 

We  now  perform  a  validation  study  analogous  to  the  DPD  fluids  study  for  a  coarse-grain  solid 
model,  where  we  consider  a  recently  developed  nickel  model  that  reasonably  reproduces  several 
measured  properties,  including  the  melting  temperature  (13).  For  a  benchmark  system,  a 
constant-temperature  DPD  simulation  is  performed  at  p  =  8260  kg/nr  and  T  =  1300  K  for 
tnm  =  1  ns  and  At  =  5  fs,  where  results  are  listed  in  table  3.  (Since  the  coarse-grain  solid  model 

has  been  parameterized  to  an  actual  material,  results  are  reported  in  real  units  as  opposed  to 
reduced  units  for  the  DPD  fluid.)  At  this  state  point,  the  atomistic  Sutton-Chen  (SC)  model  of 
nickel  predicts  a  pressure  of  approximately  0  bar  (13),  while  (Pvir)  for  the  coarse-grain  solid 

model  is  larger  than  0  bar. 

Analogous  to  the  DPD  fluid  study,  DPD-E  simulations  were  perfonned  at  conditions  determined 
from  the  constant-temperature  DPD  simulation.  The  value  of  E  imposed  in  the  DPD-E 
simulation  was  determined  from  the  final  configuration  of  the  constant-temperature  DPD 
simulation,  which  is  used  as  the  starting  configuration.  The  value  of  ui  was  initialized  by  setting 

ui  =  CvTini ,  where  Tini  -T  =  1300  K.  The  simulation  was  run  for  trun  =  1  ns  and  At  =  5  fs,  where 
a  relative  drift  in  E  under  2  •  10  4  was  observed.  (For  the  coarse-grain  solid  simulations, 
reported  relative  drifts  refer  to  an  average  of  relative  drifts  over  time  periods  of  1  ns.)  Comparing 
the  DPD-E  results  with  the  constant-temperature  DPD  results  in  table  3,  excellent  overall 
agreement  is  found.  For  DPD-E,  the  values  of  (Tkln)  and  (Tml)  are  equal  within  statistical 

uncertainties.  These  values  are  approximately  2%  lower  than  Tini  =1300  K  but  agree  within 
statistical  uncertainties  when  the  extra  degree  of  freedom  due  to  the  fluctuations  in  ui  is 
considered,  i.e.,  (Tkin)  =  (7jnt)  =  Tini  (l  -  kB  /  Cv)  =  1273  K  (16).  Furthermore,  because  of  these 
lower  values  of  (Tkin )  and  (Tm  ) ,  the  value  of  (Pvtr )  for  DPD-E  in  table  3  differs  accordingly 
from  ( Pvjr )  for  constant-temperature  DPD. 

Table  3.  The  molar  configurational  energy  (a  '"j  ,  the  kinetic  temperature 
( Tkin  ^  ,  the  internal  temperature  (Tjnt  ^ ,  and  the  virial  pressure 
(  P  jr  \  determined  from  test  case  1  simulations  of  the  coarse-grain 
solid  model  of  nickel.  (.')  denotes  an  ensemble  average,  where 

numbers  in  parentheses  are  uncertainties  calculated  from  block 
averages. 


Variant 

(u) 

(kJ/mol) 

(4.) 

(K) 

(t.  ) 

X-4  int  / 

(K) 

A) 

(bar) 

DPD 

p  =  8260  kg/m3 

-508.44(11) 

1300.1(91) 

— 

5.91(94) 

DPD-E 

p  =  8260  kg/m3 

-508.71(11) 

1274.8(88) 

1274.5(5) 

-150.49(98) 

13 


4.2  Test  Case  2:  Heating  Response  in  DPD-E  Simulations 

As  a  second  test  case  to  verify  the  SSA-VV  formulation  for  the  DPD-E  variant,  a  nonequilibrium 
scenario  was  considered  for  both  the  DPD  fluids  and  the  coarse-grain  solid  model.  Starting  from 
a  final  configuration  of  a  constant-temperature  DPD  simulation,  we  instantaneously  heated  a  slab 
of  DPD  particles  in  the  middle  of  the  simulation  box  and  studied  the  system  response  at  constant- 
(V,E)  conditions,  i.e.,  by  DPD-E  simulations. 

4.2.1  DPD  Fluid 

Analogous  to  test  case  1,  the  final  configuration  from  the  constant-temperature  DPD  simulation 
(at  T  =  1  and  p  =  3  )  was  used  as  the  starting  configuration.  For  this  configuration,  a  0.5Z  wide 

slab  of  particles  in  the  middle  of  the  simulation  box  was  heated  by  assigning  velocities  from  a 
Maxwell-Boltzmann  distribution  corresponding  to  Theat  and  by  setting  ui  =  CV  iTheat .  The 

remaining  (nonheated)  particles  were  assigned  u,  =  CV  lTini ,  where  Tjni  -  T  =  1 .  As  a  test  of  the 
DPD-E  variant,  a  simulation  was  performed  using  Theat  =  10 ,  for  tnm  -  5000  and  At  =  0.005  for 
the  pure  and  equimolar  binary  DPD  fluids.  Results  for  the  time  evolution  of  Tkin ,  Tm  and  Pvir  for 

the  pure  DPD  fluid  are  displayed  in  Fig.  1,  where  a  relative  drift  of  2  •  10  in  E  was  observed. 
From  figure  1,  it  is  evident  that  Tkjn  and  Tint  quickly  equalized,  after  which,  Tkni ,  Tm  and  Pir 

increased  with  t  until  the  system  reached  equilibrium  at  t  >  1 50 .  The  inset  of  figure  1  displays 
the  early  time  behavior  of  Tkin ,  Tml ,  and  Pvir .  Within  t  =  10  ,  Tkjn  and  Pvir  sharply  decreased  from 

their  initial  values,  followed  by  increasing  values  as  the  system  moved  toward  an  equilibrium 
state.  The  initial  dramatic  decreases  in  Tkin  and  Pvir  are  associated  with  the  relaxation  of  the 

interfaces  between  “cold”  and  “hot”  regions  in  the  simulation  box,  which  were  artificially 
created  by  the  instantaneous  heating  at  t  -  0 .  We  found  analogous  results  for  the  equimolar 
binary  DPD  fluid  (not  shown). 
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Figure  1.  Time  evolution  of  the  kinetic  temperature  Tkin,  internal  temperature  Tmt,  and  virial  pressure 
Pvir  for  a  DPD-E  simulation  of  the  pure  DPD  fluid  at  p  =  3 ,  where  a  slab  of  particles  in  the 
simulation  box  was  instantaneously  heated  by  Theat  =  1 0  at  t  =  0.  Inset  of  figure  displays 
early  time  behavior  of  Tkim  Tml ,  and  Pvir. 

4.2.2  Coarse-Grain  Solid 

A  validation  study  analogous  to  the  DPD  fluid  study  was  carried  out  for  the  coarse-grain  solid 
model  of  nickel.  The  final  configuration  from  the  constant-temperature  DPD  simulation  (at 
T  =  1300  K  and  p  =  8260  kg/m  )  was  used  as  the  starting  configuration.  From  this  starting 

configuration,  a  0.5Z  wide  slab  of  particles  in  the  middle  of  the  simulation  box  was  heated  by 
assigning  velocities  from  a  Maxwell-Boltzmann  distribution  corresponding  to  Theat  and  by 
setting  ui  =  CvTheat .  The  remaining  (nonheated)  particles  were  assigned  ui  =  CvTini ,  where 
Tini  -T  =  1300  K.  As  a  test  of  the  DPD-E  variant,  simulations  were  performed  at  p  =  8260 
kg/m3,  using  Theat  =  3000  K  for  trun  =  1  ns  and  At  -  5  fs.  We  observed  relative  drifts  of  2  ■  1 0  4 

in  E  .  At  low  and  moderate  pressures,  the  coarse-grain  solid  model  melts  between  1800  and  1850 
K  (13).  As  a  result,  at  the  end  of  the  DPD-E  run,  the  particle  configuration  corresponds  to  a 
liquid  state.  Figure  2  shows  the  time  evolution  of  Tkjn ,  Tml ,  and  Pvir  (together  with  a  few 

representative  simulation  snapshots)  for  the  DPD-E  simulation.  Complete  melting  is  evidenced 
by  reaching  a  plateau  in  the  time  evolution  of  Pvir  for  the  DPD-E  simulation,  where  complete 

melting  occurs  at  ~0.7  ns. 
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Figure  2.  (a)  Time  evolution  of  the  kinetic  temperature  Tkin,  internal  temperature  Tmh  and  virial  pressure  Pvir,  along 
with  (b)  a  few  representative  simulation  snapshots  for  a  DPD-E  simulation  of  the  coarse-grain  solid  at 
p  —  8260  kg/m3,  where  a  slab  of  particles  in  the  simulation  box  was  instantaneously  heated  by 


Them  =  3000  K  at  /  =  0. 
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Figure  2.  (a)  Time  evolution  of  the  kinetic  temperature  Tkim  internal  temperature  Tmh  and  virial  pressure  Pvir, 
along  with  (b)  a  few  representative  simulation  snapshots  for  a  DPD-E  simulation  of  the  coarse-grain 
solid  at  p  —  8260  kg/m3,  where  a  slab  of  particles  in  the  simulation  box  was  instantaneously  heated  by 
Theat  =  3000  K  at  t  =  0  (continued). 
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4.3  Conservation  of  Total  System  Energy 

In  practice,  the  SSA  involves  integrating  the  fluctuation-dissipation  contribution  and  the 
conservative  contribution  in  separate,  independent  steps.  First,  by  implementing  equation  18b 
rather  than  numerical  discretization  of  equation  5,  the  integration  of  the  fluctuation-dissipation 
contribution  exactly  conserves  the  total  energy  E  at  each  time  step.  Second,  analogous  to  an 
application  of  the  velocity-Verlet  algorithm  in  microcanonical  molecular  dynamics  (MD),  the 
integration  of  the  conservative  contribution  does  not  conserve  E  at  each  time  step.  (DPD-E 
reduces  to  microcanonical  MD  in  the  limit  of  vanishing  dissipative  and  random  forces  and  heat 
transfers.)  Rather,  the  velocity-Verlet  algorithm  preserves  E  only  up  to  terms  of  order  At2 , 
conserving  a  pseudo-FIamiltonian  that  differs  from  the  true  Flamiltonian  by  this  difference  of 
order  At2  (22-25).  Although  the  velocity-Verlet  algorithm  is  area-preserving,  it  is  not  exactly 
symplectic.  (An  algorithm  is  area-preserving  if  drdp  =  const. ,  where  r  =  jr  }^  and  p  =  {p,  }^  .) 

The  velocity-Verlet  algorithm  in  microcanonical  MD  thus  produces  a  long-tenn  energy  drift. 
Nonetheless,  because  of  its  area-preserving  property,  the  velocity-Verlet  algorithm  is  more  stable 
at  long  times  than  non- area-preserving  schemes,  since  system  trajectories  in  phase  space  that  are 
initially  close  will  remain  close  during  the  microcanonical  simulation  (22). 

For  values  of  At  comparable  to  those  used  in  constant-temperature  DPD  simulations,  we 
observed  a  small  long-term  drift  in  E  for  the  DPD-E  simulations.  For  example,  for  the  values  of 
At  used  for  the  DPD-E  simulations  in  this  work  (At  -  0.01  for  the  DPD  fluids  and  At  -5  fs  for 
the  coarse-grain  solid),  the  small  relative  drift  produced  in  E  was  typically  of  order  10  4 .  When 
At  was  decreased  by  an  order  of  magnitude,  the  relative  drift  in  E  dropped  to  order  10~7 .  A 
typical  example  of  the  dependence  of  the  relative  drift  in  E  on  At  for  the  DPD  fluid  is  shown  in 
figure  3.  The  values  of  other  properties  (not  shown  here),  such  as  the  kinetic  and  internal 
temperatures,  configurational  energy,  and  virial  pressure,  change  with  At  by  less  than  0.5%. 

This  behavior  is  comparable  with  microcanonical  MD  simulations  when  the  velocity-Verlet 
algorithm  is  used  to  integrate  the  EOM  (21).  Since  the  integration  of  the  fluctuation-dissipation 
contribution  exactly  conserves  the  energy  (up  to  machine  precision),  the  drift  is  caused  by  the 
velocity-Verlet  algorithm  during  the  integration  of  the  deterministic  contribution  in  the  DPD-E 
EOM.  Similar  to  microcanonical  MD,  a  long-term  drift  in  E  is  thus  inevitable  in  the  SSA  for 
DPD-E.  Notably,  a  recent  application  of  the  standard  velocity-Verlet  algorithm  to  DPD-E  for  the 
DPD  fluid  model  by  Abu-Nada  (7,  8)  required  At  -  0.00002  to  0.00005,  i.e.,  values  of  At  that 
are  several  orders  of  magnitude  smaller  than  for  constant-temperature  DPD  to  minimize  a  drift 
in  E . 

To  enforce  constant  energy  beyond  this  small  drift,  one  can  apply  the  following  numerical 
procedure.  After  each  time  step,  the  difference  between  the  current  E  and  the  inputted  E  is 
calculated.  This  difference  is  then  divided  by  the  number  of  particles  and  equally  subtracted  from 
each  Uj .  This  is  a  useful  strategy  provided  that  the  drift  in  E  has  a  mechanical  origin,  which 
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implies  that  the  energy  drift  scales  as  kBT  .  Thus,  the  extra  energy  per  particle  subtracted  in  this 
procedure  is  very  small  compared  to  the  magnitude  of  ut ,  which  scales  as  CT  .  In  this  work,  the 

variation  of  the  system  temperature  due  to  this  drift  was  negligible  and  the  dynamics  unaffected. 
This  strategy  was  applied  to  all  test  cases  for  DPD-E,  where  we  observed  no  variation  in  the 
results. 
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Figure  3.  The  relative  drift  in  E  as  a  function  of  the  integration  time  step  At  for 
DPD-E  simulations  with  the  SSA-VV. 


5.  Conclusion 


We  presented  a  comprehensive  description  of  a  numerical  integration  scheme  based  upon  the 
SSA  for  the  DPD-E  approach,  where  it  was  readily  extendable  and  found  to  be  a  stable  and 
accurate  integration  scheme.  The  DPD-E  variant  was  verified  using  both  a  standard  DPD  fluid 
model  and  a  coarse-grain  solid  model,  where  thermodynamic  quantities  as  well  as  probability 
distributions  were  considered.  The  integration  algorithm  for  the  DPD-E  variant  was  further 
verified  by  considering  equilibrium  and  nonequilibrium  simulation  scenarios.  Finally,  we 
discussed  the  inevitable  small,  long-term  drift  in  E  associated  with  finite  integration  methods, 
where  we  proposed  a  simple  strategy  to  minimize  the  effect  of  this  drift  in  DPD-E  simulations. 
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Implementing  the  SSA  for  a  given  conservative  force  potential,  we  found  that  a  smaller  time  step 
is  required  for  a  DPD-E  simulation,  relative  to  the  time  step  of  a  constant-temperature  DPD 
simulation.  This  behavior  is  consistent  with  the  analogs  of  MD  integration  algorithms.  Moreover, 
the  relative  sizes  of  the  time  steps  of  constant-temperature  DPD  versus  DPD-E  simulations  are 
comparable  to  the  relative  sizes  of  the  time  steps  for  canonical  versus  microcanonical  MD 
simulations.  Finally,  and  perhaps  most  importantly,  compared  to  standard  DPD  integrators  (17), 
while  the  SSA  allows  for  modest  increases  in  the  size  of  the  time  step  for  constant-temperature 
DPD  simulations,  the  SSA  allows  for  much  larger  time  steps  for  DPD-E  simulations.  Comparing 
with  the  recent  study  of  Abu-Nada  that  used  the  standard  velocity- Verlet  algorithm  for  DPD-E 
simulations  of  the  DPD  fluid  model  (18, 19),  the  SSA  proposed  in  this  work  allows  for  time 
steps  as  much  as  1 0  times  larger,  which  is  an  essential  improvement  for  practical  applications  of 
the  DPD-E  method.  So  while  the  computational  cost  of  the  SSA  is  almost  twice  that  of  the 
standard  velocity- Verlet  algorithm,  this  cost  is  compensated  by  allowing  larger  time  steps. 
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Here,  we  summarize  the  Fokker-Planck  equation  (FPE)  and  outline  the  derivation  of  the 
fluctuation-dissipation  theorem  (FDT)  for  the  dissipative  particle  dynamics  method  at  constant 
energy  (DPD-E)  variant  considered  in  the  report.  The  FPE  corresponding  to  the  equations  of 
motion  (EOMs)  given  by  equation  5  of  the  report  is 

%  =  LCP  +  LdP  +  KondP  >  (A- 1 ) 


where  the  conservative  operator  Lc  is  given  by 
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The  operator  representing  the  effects  of  the  dissipative  and  random  forces  is  given  by 
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while  the  operator  associated  with  the  effects  of  the  mesoscopic  heat  transfer  between  particles  is 
given  by 
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with  the  condition 
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=  0  to  ensure  that  the  EOMs  contain  no  spurious  drift. 


In  equations  A-l  through  A-4,  p  =  p(r,p,u;t)  and  u  =  {m.}^  are  the  particle  internal  energies; 
Ky  and  atj  are  the  mesoscopic  thennal  conductivity  and  noise  amplitude  between  particle  i  and 
particle  j ,  respectively;  coDq  and  coRq  are  weight  functions  associated  with  the  mesoscopic  heat 


transfer  between  particles;  and  6j  =  — -  is  the  internal  temperature  ( st  is  the  mesoscopic  entropy 

8s  t 


of  particle  i ). 


In  contrast  to  constant-temperature  DPD,  an  implicit  heat  reservoir  is  not  present  in  constant- 
energy  DPD.  Furthennore,  the  FPE  (A-l)  does  not  impose  external  constraints  on  the  system 
(such  as  the  total  system  energy).  Thus,  for  the  purposes  of  deriving  the  FDT,  it  is  equivalent  to 
consider  that  the  system  is  either  isolated  or  in  thermal  contact  with  a  heat  reservoir,  while  the 
resulting  FDT  relations  for  either  choice  should  be  insensitive  to  any  external  constraints.  The 
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use  of  the  canonical  distribution  simplifies  this  derivation.  Therefore,  if  we  consider  that  the 
system  is  in  contact  with  a  heat  reservoir  that  maintains  the  system  temperature  T ,  then 
peq  =  peq{r, p,  u)  corresponds  to  the  canonical  probability  density 
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(A- 5) 


Similarly  as  before,  Lcpeq  =  0 , 1  while  the  FDT  then  follows  from  the  requirements  that 
LDpeq  =  0  and  Lcondpeq  =  0 .  LDpeq  -  0  is  satisfied  for 

=2rij(oDkB&i], 


i.e.,  for 
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while  Lcondpeq  =  0  is  satisfied  for 


i.e.,  for 


al  = 2  h*ij 


(A-6) 


(A-7) 


Generally,  /c:j  can  be  a  function  of  u,  and  u  - .  For  such  a  case,  Lcondpeq  =  0  requires  that 
Ky  =  /c[uj  +  Uj  )  because  of  the  constraint  on  a*  given  in  equation  A-4.  As  expected,  the  FDT 

relations,  equations  A-6  and  A-7,  do  not  depend  on  the  heat  reservoir  temperature  and  therefore 
do  not  depend  on  the  ensemble  used  for  the  derivation  of  these  relations. 


1  Brennan,  J.  K.;  Lisal,  M.  Dissipative  Particle  Dynamics  at  Isothermal  Conditions  Using  Shardlow-Like  Splitting  Algorithms', 
ARL-TR-6582;  U.S.  Anny  Research  Laboratory:  Aberdeen  Proving  Ground,  MD,  September  2013. 
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Intentionally  left  blank. 
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For  the  models  considered  in  this  work,  the  details  of  the  conservative  forces  expressed  in 
equation  2  of  the  main  text  are  the  following.  u^G  for  the  pure  and  binary  Dissipative  Particle 


Dynamics  (DPD)  fluid  models  is  given  by 

UT  =  avrca>D(rv), 


(B-l) 


where  atj  is  the  maximum  repulsion  between  particle  i  and  particle  j . 


For  the  coarse-grain  solid  model,  which  has  a  face-centered-cubic  (f.c.c.)  lattice  structure, 
particles  interact  through  a  shifted-force  Sutton-Chen  embedded  potential  (SC)  embedded 
potential  given  as: 


CG 

U,  =  £ 


1 


V  Z  J*i  J 


'4f>i 


(B-2) 


where 


(B-3) 


e  and  r0  are  the  energy  and  length  parameters,  respectively,  n  and  m  are  positive  integers 
(n>  m  to  satisfy  elastic  stability  of  the  crystal),  and  c  is  a  dimensionless  parameter.  Although 
effectively  this  is  a  many-body  potential,  the  force  on  each  particle  can  be  written  as  a  sum  of 
pairwise  contributions.  The  coarse-grain  solid  model  used  here  approximates  nickel,  where  one 
DPD  particle  was  chosen  to  represent  four  f.c.c.  unit  cells,  i.e.,  16  nickel  (Ni)  atoms.  SC  potential 
parameters  were  determined  by  fitting  to  various  0-K  properties  and  the  melting  temperature  at 
zero  pressure,1 2 3  where  the  following  values  were  found:  £t  kB  =  225  K,  r0  -  8.8698  A, 

c  -  39.43 14,  m  =  6 ,  and  n  =  9  .  Further  details  for  determining  SC  parameters  based  upon  such 
a  procedure  can  be  found  elsewhere.2,3 


1  WebElements.  http://www.webelements.com/  (accessed  16  August  2013). 

2 

“Brennan,  J.  K.;  Llsal,  M.  Proceedings  of  the  14th  International  Detonation  Symposium ,  Coeur  d’Alene,  ID,  1 1-16  April 
2010. 

3Sutton,  A.  P.;  Chen,  J.  Long-Range  Finnis  Sinclair  Potentials.  J.  Philos.  Mag.  Lett.  1990,  61,  139. 
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List  of  Symbols,  Abbreviations,  and  Acronyms 


DPD 

constant-temperature  Dissipative  Particle  Dynamics 

DPD-E 

constant-energy  Dissipative  Particle  Dynamics 

EOMs 

equations  of  motion 

f.c.c. 

face-centered-cubic 

FDT 

fluctuation-dissipation  theorem 

FPE 

Fokker-Planck  equation 

MD 

molecular  dynamics 

SC 

Sutton-Chen  embedded  potential 

SDEs 

stochastic  differential  equations 

SSA 

Shardlow-splitting  algorithm 

SSA-VV 

Shardlow-splitting  algorithm-velocity  Verlet 

29 


NO.  OF 

COPIES  ORGANIZATION 

1  DEFENSE  TECHNICAL 
(PDF)  INFORMATION  CTR 
DTIC  OCA 

1  DIRECTOR 

(PDF)  US  ARMY  RESEARCH  LAB 
IMAL  HRA 

1  DIRECTOR 

(PDF)  US  ARMY  RESEARCH  LAB 
RDRL  CIO  LL 

1  GOVT  PRINTG  OFC 
(PDF)  AMALHOTRA 

ABERDEEN  PROVING  GROUND 

1  DIR  USARL 
(PDF)  RDRL  WML  B 
J  BRENNAN 


30 


